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Effects of thermal and quantum fluctuations on the phase diagram of a spin-1 87 Rb 
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We investigate effects of thermal and quantum fluctuations on the phase diagram of a spin-1 87 Rb 
Bose-Einstein condensate (BEC) under a quadratic Zeeman effect. Due to the large ratio of spin- 
independent to spin-dependent interactions of 87 Rb atoms, the effect of noncondensed atoms on the 
condensate is much more significant than that in scalar BECs. We find that the condensate and 
spontaneous magnetization emerge at different temperatures when the ground state is in the broken- 
axisymmetry phase. In this phase, a magnetized condensate induces spin coherence of noncondensed 
atoms in different magnetic sublevels, resulting in temperature-dependent magnetization of the 
noncondensate. We also examine the effect of quantum fluctuations on the order parameter at 
absolute zero, and find that the ground-state phase diagram is significantly altered by quantum 
depletion. 

PACS numbers: 03.75.Hh,03.75.Mn,67.85. Jk 



I. INTRODUCTION 

Since the first experimental realization of a Bose- 
Einstein condensate with spin degrees of freedom (spinor 
BEC) in 1998 [HQ, many interesting phenomena have 
been investigated. Due to the competition between the 
interatomic interactions and the coupling of atoms to 
an external magnetic field d, H| , these systems can ex- 
hibit various phases having different spinor order pa- 
rameters 0. Both theoretical and experimental stud- 
ies have extensively been conducted on various aspects 
of spinor BECs (see, for example, Q). Experiments 
have been performed to investigate formation of spin do- 
mains Q or tunneling between them 

0- Spin -mixing 

dynamics has also been observed in both spin-1 and 
spin-2 BECs [8T-[l2j. More recently, precise control of 
the magnetic field has enabled experimenters to ob- 
serve amplification of spin fluctuations (l3l . [h"H and real- 
time dynamics of spin vortices and short-range spin tex- 
tures |15l - ll7| . Finite-temperature properties of spinor 
BECs have also been theoretically investigated: the dy- 
namics of spinor systems in quasi-one [H|,ll9| and three- 
dimensional spaces [20| , and finite-temperature phase di- 
agrams of both ferromagnetic |2l|-|23j and antiferromag- 
netic spinor condensates [2J-|26 

For scalar BECs, the first-order self-consistent approxi- 
mation (also called the Popov approximation [27j ) , which 
neglects the pair correlation of noncondensed atoms or 
the anomalous average, can give a good description of 
thermal equilibrium properties of the system over a wide 
range of temperatures except near the BEC transition 
point. This is because at temperatures well above abso- 
lute zero, the anomalous average is negligibly small com- 
pared with the noncondensate number density. In con- 
trast, near absolute zero the anomalous average is of the 
same order of magnitude as the noncondensate number 
density, but both of them are very small compared with 



the condensate density and hence negligible. However, 
for spinor BECs, in particular, spin-1 87 Rb BECs, due 
to the large ratio of spin- independent to spin-dependent 
interactions, the anomalous average and noncondensate 
number density are expected to be as important as the 
spin-dependent interaction between two condensed atoms 
near absolute zero. 

The above striking difference between the scalar and 
spinor BECs has hitherto not been fully studied. A full 
investigation of this problem is the main theme of this 
paper. In Refs. (2lL |2^|. the quadratic Zeeman energy, 
which is a key control parameter in spinor BECs, was not 
taken into account. In the present theoretical study, we 
investigate effects of thermal and quantum fluctuations in 
a spinor Bose gas in the presence of the quadratic Zeeman 
effect. We consider a three-dimensional uniform system 
of spin-1 87 Rb atoms with a ferromagnetic interaction, 
where the spin-independent interaction is stronger than 
the spin-dependent interaction by a factor of about 200. 
Therefore, even when the fraction of noncondensed atoms 
is small, they can significantly affect the magnetism of the 
system via the spin-independent interaction. 

In this paper, we first use the first-order self-consistent 
approximation to obtain the finite-temperature phase di- 
agram in the presence of a quadratic Zeeman effect. We 
find that the system undergoes a two-step phase tran- 
sition, where condensation and spontaneous magnetiza- 
tion occur at different temperatures. We then examine 
temperature-dependent magnetization of the nonconden- 
sate, which is a remarkable consequence of the spin co- 
herence induced by the magnetized condensate. To in- 
vestigate the effect of quantum depletion on the phase di- 
agram at absolute zero, we adopt the method developed 
in Ref. [28[ , in which the order parameter is expanded in 
powers of the square root of the noncondensate fraction. 
By applying the method to spinor systems, we find a sig- 
nificant modification of the ground-state phase diagram 
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due to the effect of noncondcnscd atoms. 

This paper is organized as follows. Section [TT] intro- 
duces a theoretical framework of spin-1 spinor BECs, and 
describes the mean-field ground-state phase diagram an- 
alytically. Section IIIII discusses the finite-temperature 
phase diagram by using the first-order self-consistent ap- 
proximation, and studies magnetizations of the conden- 
sate and noncondensate as functions of temperature. Sec- 
tion IIVI investigates the effect of quantum depiction on 
the zero-temperature ground-state phase diagram. The 
perturbative expansion method for spinor BECs is intro- 
duced, followed by a discussion of a modification of the 
ground-state phase diagram from the first-order counter- 
part. Finally, Sec. fVl concludes this paper by discussing 
possible experimental situations. Complicated algebraic 
manipulations that would distract readers from the main 
subject arc placed in Appendices. 



poncnts of the spin-1 matrix vector given by 



II. HAMILTONIAN AND MEAN-FIELD 
GROUND STATE 



We consider a system of spin-1 identical bosons with 
mass M that are confined by an external potential U (r) 
and subject to a magnetic field in the z direction. The 
one-body part of the Hamiltonian is given in matrix form 
by 



(ho)i 



2V72 



h 2 \7 

2M 



+ U(r) — pi + qv 



(1) 



where the subscripts i,j = 0, ±1 refer to the magnetic 
sublevels, and p and q are the coefficients of the linear and 
quadratic Zccman terms, respectively. The total Hamil- 
tonian of the spin-1 spinor Bosc gas is given in the second 
quantization by 0, \ 



H 



$( r )(My4-( r ) 



-^^t( r )^t(r)4-(r)^i(r) 



+ y £ (/a)«(/a)w^ t (r)^(r)^(r)^(r) > (2) 



where ipi (r) is the field operator that annihilates an atom 
in the magnetic sublevel i at position r, a = x, y, or z 
specifies the spin component, and / Q 's denote the com- 
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(3) 
(4) 
(5) 



The last two terms in the Hamiltonian (J2J describe the 
spin-independent and spin-dependent interactions, re- 
spectively. The coefficients Co and c\ can be expressed 
in terms of the s-wave scattering lengths ao and ci2 of 
binary collisions with total spin F to t a i = and 2, respec- 
tively, as [1] 



Co 

Cl 



Anh 2 a + 2d2 
:_ M 3 : 

Airti 2 ci2 — ao 
'~M 3 ' 



(6a) 
(6b) 



In the mean-field ground state of a spinor Bose gas, the 
effect of quantum depletion is neglected, and all particles 
are assumed to occupy the same single-particle state in 
both coordinate and spin spaces. The field operator ^i(r) 
can then be replaced by a classical field <fri(r), and the 
expectation value of Hamiltonian ([2]) is given by the fol- 
lowing energy functional: 



EM = dr 



(7) 



where the number density n c (r) and the three compo- 
nents of the spin density vector F c (r) of the condensate 
are given by 

n^^l^r)! 2 , (8) 

i 

JS(r)=£#(r)(/ a )«fc(r) (a = x,y,z). (9) 

In the mean-field approximation, n c is equal to the total 
number density n. The condensate wave function 4>i(r) 
is determined by minimizing the energy functional (J7]), 
i.e., 



5E[<j>i 



0. 



subject to the normalization condition 



J drJ2\Mr)\ 2 = N, 



(10) 



(11) 
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where N is the total number of atoms. Equation (TT01) . 
together with Eq. (fTTj) . leads to the multi-component 
Gross-Pitaevskii (GP) equation: 



E 



(ho) 



c n8i 



ci£^(/ Q 



<t>3=lulH, (12) 



where fi is the chemical potential at absolute zero. 

For a uniform system, i.e., when U(r) = 0, the con- 
densate wave function is independent of r and the 
solutions to Eq. (|T2l) can be obtained analytically. For 
the case of c\ < and p = 0, which is the case we con- 
sider in the present paper, the order parameters = 
(<j>i, <f)Q, 4>-i) T and the energies per particle e = E[<j>j\/N 
for possible phases are given as follows 0, [29| : 



Ferro : = V^(1,0,0) T or Vn(0, 0, 1) T , 

CO + Cl 



e = 5 



-n, 



Polar : cj> = V^(0,1,0) T , 
Co 



/ 



BA : 4> = \ — 



1 



-iO 



^0 



2\ Cl \n 



2\ci\n 



l - 



2 cxln 



Cl 



2|ci|n 

c 



-n, 



(13) 
(14) 
(15) 
(16) 

(17) 
(18) 



where T denotes transpose, and Ferro, Polar, and BA 
stand for ferromagnetic, polar, and broken-axisymmetry 
phases, respectively. In Eq. (fl~7|) . 8 can take on values 
between and 2w, and we have omitted overall phase 
factors in Eqs. (fl3]l. (fl~5j) . and ([171). The BA phase exists 
only in the region of < q < 2\ci\n, and becomes the 
ground state of the system in this parameter regime. The 
magnetization for the BA phase is transverse and given 
by 



F c = 



r + — r X 



iF* = ne' 



1- 



2cin 



(19) 
(20) 



Hence, 9 specifies the direction of magnetization in the 
xy-plane, and its magnitude depends on q. The BA phase 
is named after the fact that the transverse magnetiza- 
tion breaks the rotational symmetry of the Hamiltonian 
around the z-axis [2j|. If q > 2|ci|n, the ground state 
is in the polar phase. In this phase, the condensate has 
zero magnetization. On the other hand, if q < 0, the 
fully polarized state in the magnetic sublevel i = 1 or 
— 1 minimizes both the ferromagnetic interaction and the 
quadratic Zecman energy. Therefore, the ferromagnetic 
phase is the ground state of the system. To satisfy the 
conservation of the total longitudinal magnetization, a 




FIG. 1: (Color online) q dependence of longitudinal mag- 
netization per particle F^/n c and that of transverse one 
F1/n c in the mean-field ground state of a spin-1 ferromag- 
netic BEC (ci < 0). The dashed line and solid curve show 
the longitudinal and transverse magnetizations, respectively, 
where F| = = s/(Fg)* + (F°) 2 . Ferro, BA and Po- 

lar stand for ferromagnetic, broken-axisymmetry, and polar 
phases, which are shaded in light, medium and dark grey, 
respectively. The longitudinal magnetization per particle 
is F^/n c — Q(—q), and the transverse one is |F£|/n c = 
yjl - (q/2|ci|n) 2 e(<?)e(2|ci|n - q), where 0(q) is the unit- 
step function. 



phase separation into two spin domains with F^/n c = 1 
and —1 must occur. Figure [T] shows the q dependence of 
the longitudinal magnetization and that of the transverse 
one. 



III. FINITE-TEMPERATURE PHASE 
DIAGRAM UNDER THE FIRST-ORDER 
SELF-CONSISTENT APPROXIMATION 

A. First-order self-consistent approximation 

At finite temperatures, a fraction of atoms are ther- 
mally excited from the condensate to form a thermal 
cloud, which, in turn, will affect the condensate. There- 
fore, the finite-temperature phase diagram should be de- 
termined in a self-consistent manner. The field operator 
is decomposed into the condensate part, which can be 
replaced by a classical field 4>i(r), and the noncondensate 
part Si(r): 



(21) 



For convenience, we consider here a grand-canonical en- 
semble of the atomic system and introduce the operator 



K, = U-iMAf, 



(22) 
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where the total number operator JV is denned as 



d* £#(')&(*)• 



(23) 



Substituting Eq. $ZT§ into Eqs. © and ([23]), and collect- 
ing terms of the same order with respect to the fluctua- 
tion operator Si(r), we obtain 



fC = K + K 1 + K 2 + K : - 



(24) 



where K n (n = 0, . . . , 4) is comprised of the terms that 
involve the n-th power of Si(r). 

The static properties of the system in thermal equi- 
librium can be calculated from the eigenspectrum of op- 
erator /C. The part of K, that involves the terms up to 
quadratic in Si, i.e., K0+K1+K2, can be diagonalized us- 
ing a Bogoliubov transformation [30j . The higher-order 
terms in K3 and K4 can be made into quadratic forms by 
applying the mean-field approximation to noncondensed 
atoms. The mean-field approximation for noncondensate 
operators in K 3 and if 4 is carried out as follows [27l . [3l| : 



SjSjSk ~ hijSk + n t k5j + 6[rhjk, 
Si ~ ^45; + nji5\5k - 
Si 

+ rhijSkSi + fhuS\5] - rh^rhki, 



+ nuS]S k + n jk 5\Si - huh jk 



(25) 



(26) 



where JVj(r) = (i5](r)<5j(r)) is the matrix element 
of the noncondensate number density, and rhij (r) = 
(Si(r)Sj(r)) is that of the noncondensate pair correlation, 
which is called the anomalous average. 

In the first-order self-consistent approximation, the 
anomalous averages rhij ( r ) an d rh*j (r) are neglected [27| ■ 
This would give a gapless spectrum of elementary excita- 
tions, in agreement with the Nambu-Goldstone theorem. 
At temperatures well above absolute zero, the spectrum 
of elementary excitations approaches that of single parti- 
cles, and, therefore, the anomalous average rhij becomes 
negligibly small compared with the noncondensate num. 



where the condensate number density n° (r) and spin den- 
sity -F^(r) (a = x 7 y,z) are defined in Eqs. © and ©, 
respectively, and n nc (r) and F£ c (r) are the nonconden- 
sate counterparts given by 



n nc (r)^]TMr), 

i 

F^r)^(/ Q )^(r). 



(28) 
(29) 



Thus, the operator K, reduces to the sum of a c-number 
termi^o and a quadratic operator K.^ , where 



Kn 



dr 



|(n c ) 2 + ^|F c | 2 



5\Aij{v)5j + 



(30) 




\\5\Bij{v)5] 



(31) 



Here, the matrices Aij(r) and -B,j(r) are defined as 
Aij{v) = (h )ij - fxSij + c [(?7 C + n DC )Sij 



ci 



£ 



+ ^2{fa)il(fa)kj{(p*k(t>l + hkl) 
k,l 

Bij(r) EE C (/>i4>j +Ci ^ (fa)ik(fa)jl</>k<fil- 



(32a) 
(32b) 



a,k,l 



We diagonalize the quadratic operator KS 2 ' by a Bogoli- 
ubov transformation 



"^6"6 lul J umlul ^U"i r u.iv.u vviun Mil. iiuii^wim^n.ju.u^ mini „ 

bcr density n SJ [32(. Consequently, the first-order self- U x) = / dr^] Lf ) *(r)<5 l (r) - w 4 (A) (r)(5j(r) , (33) 

ennsistent armrrndmatinn enves a mod Heserirttinn nf a J 



consistent approximation gives a good description of a 
Bose gas in thermal equilibrium over a broad range of 
temperatures, except near absolute zero. 

The condensate wave function </>i(r) then satisfies the 
generalized GP equation, which is obtained from the re- 
quirement that the operator /C be stationary with respect 
to 4>i (r) , or equivalently, that the sum of terms that are 
linear in Si(r) vanish: 

£ \ (ho)ij<t>j + co(n c + n nc )Sij + coh^^j 
j I 

(F c a + F^)(f a ). lJ + Y,(fMf a ; 



where the coefficients u[ X \r) and v^ X \r) (i = 0, ±1) sat- 
isfy the generalized Bogoliubov-de Gennes (BdG) equa- 
tion for the excitation mode labeled by index A: 



-BUv) -A* ? (r; 



A 3 {v) Bij(r) \ (uf(r)\_ w (u^>(r) 



a 



lj n kl 



(27) 



(34) 

In thermal equilibrium, the noncondensate number den- 
sity is expressed in terms of u^(r) and v^ x \r) as 

n lJ (r)=^{ U f ) *(r)4 A) (r)/( e W) 

A 

+ «f ) (r>f >)[/(e< A )) + l] }, (35) 
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where /(e) = l/[cxp(e/fcsr) — 1] is the Bose-Einstein dis- 
tribution function, and the coefficients u^(r) and i/ A )(r) 
are normalized as 



dr£[|^ A) (r)| 2 -K (A) (r) 



1. 



(36) 



Finally, the condensate and noncondensate satisfy the 
following number equation: 



N = / dr K(r) +n llc (r)] 



(37) 



B. Finite-temperature phase diagram 

In the following sections, we consider a three- 
dimensional uniform system of spin-1 87 Rb atoms with 
a fixed total number density n. Then, the condensate 
wave function ^ and the normal density fry are constant, 
while the coefficients of the Bogoliubov transformation 
are given by 



Uj (r) = e 

(A) i \ (i/,k) i 

(r) = - e 



(38a) 
(38b) 



where k is the wave vector and v is an index to distinguish 
between excitation modes. 

We consider the case in which the system is initially 
prepared so that the total magnetization projected along 
the z-axis vanishes. Due to the conservation of the total 
longitudinal magnetization, the linear Zeeman term van- 
ishes, and, therefore, we havep = 0, q ^ in Eq. (P). The 
s-wave scattering lengths of the 87 Rb atom in the F = 1 
hyperfine manifold are calculated to be ao = 101.8 as 
and a,2 — 100.4 ae [33| . where as is the Bohr radius. 
Consequently, ci given in Eq. (J6j) is negative, i.e., the 
interaction is ferromagnetic, and it is about 200 times 
smaller than cq- 

We have numerically solved a set of coupled equa- 
tions in the first-order self-consistent approximation 
[Eqs. ([27 p — (|37 p ] at a given temperature and for a given 
value of q. Here, the generalized GP equation l|27p was 
solved numerically by using the imaginary-time propa- 
gation method, which evolves a randomly chosen initial 
state to a local minimum of the Hamiltonian. 

Figure [2] shows the finite-temperature phase diagram 
of a spin-1 87 Rb BEC with (a) n = 1.0 X 10 12 cm~ 3 
and (b) n = 1.0 X 10 13 cm~ 3 . Here, the phase of the 
system is identified by calculating the condensate number 
density n c and the longitudinal F£ and transverse F]_ = 

\J (Fx) 2 + {Fy ) 2 magnetizations of the condensate. The 
high-temperature normal phase has n c — 0, while the 
condensed phases have n c ^ 0. The ferromagnetic, BA, 
and polar phases are characterized by F^/n c = 1, F]_ ^ 
0, and F% = F]_ = 0, respectively. 

Figure [2] shows that the region of BA phase shrinks 
with increasing temperature because thermal fluctua- 
tions suppress the transverse magnetization. The phase 
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FIG. 2: (Color online) Finite-temperature phase diagram of a 
spin-1 ferromagnetic Bose gas in the first-order self-consistent 
approximation. We have used the interaction parameters for 
87 Rb atoms, i.e., ao = 101.8 ae and ei2 = 100.4 cib HH, with 
the total number density (a) n = 1.0 x 10 12 cm~ 3 , and (b) 
n = 1.0 x 10 13 cm~ 3 . The quadratic Zeeman energy q and 
temperature T are measured in units of the spin-dependent 
interaction \ci\n and the transition temperature of a uni- 
form ideal scalar BEC, T = 3.31fi 2 n 2/3 /(fc B M), respectively. 
Crosses show temperature T c i below which the condensate 
density n c becomes nonzero. Open triangles show tempera- 
ture T C 2 below which the condensate acquires a nonzero trans- 
verse magnetization F]_ . The solid curves show guides for the 
eye. 



boundary between the BA and polar phases is almost in- 
dependent of the total number density n except near ab- 
solute zero. We note that if the ground state is in the BA 
phase, the phase transition is a two-step process: first, 
the system undergoes a phase transition from the normal 
phase to the polar phase at temperature T c i ~ 0.48 To, 
where To = 3.31?i 2 n 2 / 3 / (feM) is the transition tempera- 
ture of a uniform ideal scalar BEC with the same atomic 
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density n. Here, T c i is smaller than To by a factor of 
about (1/3) 2 / 3 ~ 0.48, because above T cl , the popula- 
tion of atoms in each magnetic sublevcl is almost equal to 
?i/3. The quadratic Zeeman effect causes a further small 
shift of T c i from the value 0.48 To by making a slight 
population imbalance. A slope of the normal-condensate 
phase boundary caused by the quadratic Zeeman effect 
is too small (~ 10 -4 ) to be seen in Fig. [21 At a lower 
temperature T c2 (T c2 < T c i), the system undergoes a sec- 
ond phase transition to the BA phase having a nonzero 
transverse magnetization. 

The physics of the two-step phase transition can be 
understood as follows. Due to the positive quadratic Zee- 
man energy, the magnetic sublevel i = has a higher pop- 
ulation than the levels i = ±1. Consequently, the system 
first condenses into the polar phase. If the temperature 
is further lowered, the other states (i — ±1) also undergo 
Bose-Einstein condensation, and the system enters the 
BA phase by developing transverse magnetization. 

We note that the phase diagram shown in Fig. [5] does 
not agree with the mean-field phase diagram, in which 
the phase boundary between the BA and polar phases at 
T = is given by q = % with q\, = 2\ci\n. In the first- 
order self-consistent approximation, the phase boundary 
shifts to % = 2.12|ci|n for re = 1.0 X 10 12 cm~ 3 and 
<7b = 2.37|ci|ti for n = 1.0 x 10 13 enr 3 due to quan- 
tum fluctuations. The results suggest that a more careful 
treatment needs to be made for the anomalous average 
near absolute zero. We shall discuss this point in Scc. lIVI 

C. Condensate fraction and magnetization 

Next, we study the temperature dependence of the con- 
densate fraction n c /n, and that of the longitudinal and 
transverse magnetizations per particle of both the con- 
densate F^ ± /n c and noncondensate F^ c ± /n nc . Figure [3] 
shows the result of numerical calculation for q — \ci\n. 
For other values of q in the region < q < qb, these 
physical quantities depend on temperature in a manner 
qualitatively similar to the case of q = \ci\n. The longi- 
tudinal magnetizations F^' nc vanish over the whole tem- 
perature region. The transverse magnetizations of the 

condensate FJ_ = (F£) 2 + (F£) 2 and noncondensate 
= aJ(F™) 2 + (F£ c ) 2 are given by 

F c x + iF c y = \/2(<^ + 0*<M), (39a) 

F™ + iF™ = V2(ni, + no,-i), (39b) 

where the matrix elements of the noncondensate number 
density matrix are defined as jiy = (SjSj) below Eq. (EU). 
The transverse magnetization of the noncondensate is 
parallel or antiparallel to that of the condensate, and 
we set a = 1 (a = —1) if they are parallel (antiparallel). 

Figure [3] demonstrates a two-step phase transition, in 
which a nonzero condensate fraction emerges at temper- 
ature T c i ~ 0.48 Tq, and finite transverse magnetizations 




FIG. 3: (Color online) Temperature dependence of the con- 
densate fraction n c /n (squares), transverse magnetizations 
per particle of the condensate FJ_/n c (circles) and nonconden- 
sate Fl c /n uc (triangles) for q= \ci\n and n = 1.0 xlO 12 cm" 3 . 
The BA, polar, and normal phases are shaded in medium, 
dark, and light grey, respectively. The inset shows enlarged 
behaviors of these physical quantities near absolute zero. The 
longitudinal magnetizations of the condensate and noncon- 
densate vanish for q > 0; thus, for the BA phase their mag- 
netizations are both perpendicular to the external magnetic 
field. The negative values of the transverse magnetization of 
the noncondensate at 0.01 To < T < 0.3 To imply that the 
magnetization vectors of the condensate and noncondensate 
are antiparallel to each other. 



of both the condensate and noncondensate emerge at a 
lower temperature T c2 ~ 0.3 To. The nonzero transverse 
magnetization of the noncondensate is a consequence of 
the spin coherence of noncondensed atoms, which results 
from their coupling to the magnetized condensate. Above 
T c i, where there is no condensate, no spin coherence ex- 
ists within the thermal cloud. In contrast, above T c i, 
the condensate induces spin coherence of noncondensed 
atoms as indicated by hij 7^ for i ^ j, leading to a 
nonzero transverse magnetization F^_ c . The spin coher- 
ence was experimentally observed in a two-level spinor 
system at low temperatures [34l . [35[ . 

It can also be seen from Fig. [3] that the magnetization 
of the condensate and that of the noncondensate are an- 
tiparallel to each other over a wide range of temperatures 
(0.01 < T/To < 0.3) except near absolute zero, where 
they become parallel to each other. This phenomenon 
can be understood by considering the energy spectra of 
the excitation modes of the system described by Eq. (|34p . 
From Eqs. (|3"5| and (|39b[) . the transverse magnetization 
of the noncondensate can be expressed in terms of the 
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excitation modes as 



(a) 



r + — r x 



iF™ 



EEK„ k /(^ k) ) + ^, 1 



(40) 



v k 



where v = 1,2,3 denote the index of the excitation modes 
[see Eq. {35}], and 



F 



tf 

+ ,iAk 



= V2 



fi/,k)* (i/,k) . (^,k)* (i/,k) 



(z/,k) (iAk)* (t'jk) (iAk)* 



V V, 



(I 



v 'VI 



F 



qd 



V2 



(V.k) (V,k)* , (f,k) (V,k)* 



(41a) 
(41b) 



give the contribution to F™° from the thermally excited 
collective modes and that from the quantum depletion at 
absolute zero, respectively. 

Figures SJa) and|3Jb) show the energy spectra e'" ,k ) of 
the three modes [y = 1, 2, 3) at T = and 0.2 To, respec- 
tively. Figures IDJc) and|Ud) plot F± uk 



F 



qd 



a\F 



qd 



a|F^ k | and 



-1) 



J-.i/.k - "l^+^kb respectively, where a = 1 (a 
if the magnetization of the condensate and that of the 
noncondensate are parallel (antiparallel) to each other. 
It can be seen from Fig. fjjc) that the v ~ 1 mode has 
no magnetization (dotted line), while the other two have 
magnetizations parallel (y = 2, solid curve) and antipar- 
allel [y = 3, dashed curve) to that of the condensate. 
Note here that the excitation energy of the v = 2 mode 
is higher than that of the v = 3 mode at high momenta 
[Figs. S{a) andQJb)]. Consequently, at high tempera- 
tures, the number of thermally excited quasiparticles in 
the v — 2 mode is smaller than that in the v — 3 mode, 
leading to the negative T" , which implies that the non- 
condensate is magnetized in the direction antiparallel to 
that of the condensate. The above difference between the 
energy of the v = 2 and v = 3 modes can be explained as 
follows. If the noncondensed atoms have spin configura- 
tions differing from that of the condensate (y = 3), they 
interact with the condensate only via the direct (Hartree) 
term. In contrast, if they have the same spin configura- 
tion as the condensate (y = 2), both the direct (Hartree) 
and exchange (Fock) exist, making the excitation energy 
of the v = 2 mode higher than that of the v = 3 mode. 

On the other hand, in the low-momentum regime, the 
v = 2 mode has a gapless linear dispersion relation, which 
results in a nonzero F^ d 2 k [Fig.[4|d)], whereas the v = 3 
mode has an energy gap, which suppresses the quantum 
depletion at absolute zero. Consequently, the magneti- 
zation of the noncondensate becomes parallel to that of 
the condensate in this very low-temperature region. The 
temperature at which the magnetization of the noncon- 
densate changes its direction is T q( j — 0.01 T). Below it, 
the effect of quantum depletion becomes significant. This 
crossover temperature is, however, much higher than the 
energy gap of the v = 3 mode (~ \c\\n/kB ~ 9x 10~ 5 T ), 
which is a spinor manifestation of the Bose enhancement 
in the presence of a magnetized condensate. 



0.002 




0.001 ■ 



0.005 0.01 



(b) 



£ (vM) /(k B T Q ) 
0.06 : 

r=o.2r 




(d) f<ld 




7=0 



0.5 



1 



k/k 



FIG. 4: (Color online) Energy spectra of excitation modes at 
(a) T = and (b) T = 0.2 To, and the contributions of these 
modes to the transverse magnetization T" c of the noncon- 
densate due to (c) thermal fluctuations Tl f „ k (T = 0.2 To) 
and (d) quantum depletion F^_ v k (T = 0) [see Eqs. (|40[1 and 



(Tflj) ] for q — \ci\n and n = 1 x 10 cm" . The magnitude of 
wavevector k = |kj is measured in units of ko = v^MfceTo/ft. 
There are in total three excitation modes labeled by v = 1, 2, 
and 3, which are shown by the dotted, solid, and dashed 
curves. The energy spectra of v = 1 and v — 3 modes al- 
most coincide in the high- momentum region as shown in (b). 
The negative values of F± jJ/=3 k in (c) imply that the trans- 
verse magnetization of this mode is antiparallel to that of the 
condensate. Note that T^ „ =1 k in (c) and i 7 J qd ly=1 3 k in (d) 
vanish. 
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IV. EFFECTS OF NONCONDENSED ATOMS 
ON THE GROUND STATE AT ABSOLUTE ZERO 



In dilute, weakly interacting Bosc gases, the fraction of 
noncondensed atoms due to quantum depletion at abso- 
lute zero is very small. For example, for a uniform scalar 
BEC of 87 Rb atoms with atomic density n = 10 12 cm -3 , 
the quantum depletion at absolute zero is evaluated to 
be n nc /n = 8(?ia 3 ) 1 / 2 /3 ~ 5 x 10" 4 (see [H, p. 233) 
and its effect on the condensate is only to shift the chem- 
ical potential. Even for a trapped system, such a small 
fraction of noncondensed atoms hardly affects the shape 
of the condensate. For a spinor BEC, however, the quan- 
tum depletion significantly alters the magnetism of the 
condensate as we have discussed at the end of Sec. IIII Bl 
the phase boundary between the BA and polar phases 
shifts due to quantum depletion. 

The reason why a minute quantum depletion leads to 
a significant change in the ground-state magnetism can 
be understood from the generalized GP equation (f27|). 
In the mean-field approximation, the order parameter 
of a uniform system is determined by the competition 
between the quadratic Zeeman energy ~ qi 2 Sij and the 
spin-dependent interaction ~ ciF£(f a )ij- In the first- 
order self-consistent approximation, since |ci| -C Co, the 
noncondensed atoms affect the ground-state wave func- 
tion mainly via the terms con nc and Con*j in the first 
line of Eq. (|27|) . Here, the term con Ilc merely shifts the 
chemical potential as in the case of a scalar BEC, whereas 
the term Coh*j mixes the spinor components V'j's, thereby 
changing the magnetism of the condensate. Note that for 
the 87 Rb atoms in F = 1 hyperfine manifold, the spin- 
independent interaction is about 200 times larger than 
the spin-dependent interaction. Due to such a large ratio 
co/|ci|, the term Coh*j can have a magnitude compara- 
ble to that of the spin-dependent interaction ciF^(f a )ij 
between condensed atoms even at absolute zero. 

We therefore need to investigate the effect of the 
anomalous average, which is neglected in the first-order 
self-consistent approximation, on the ground-state mag- 
netism. To take into account the effects of both the 
anomalous average and the noncondensate density, we 
use the perturbative expansion in powers of x 1 ^ where 
X = n nc /n is the noncondensate fraction which is small 
at absolute zero. This approach was first proposed by 
Castin and Dum for scalar BECs [Hf , and we generalize 
it to spinor gases. 



A/ 2 



perturbative expansion 



According to Penrose and Onsager |37[ , the condensate 
wave function, which plays the role of the order param- 
eter, is defined as the eigenfunction of the one-particle 
reduced density matrix p\Ar,r' ,t) = (ipUr')ipi(r)}t with 



a macroscopic eigenvalue: 

f dr' 5>^(r,r',i)^(r',i) = N c (t)ip t (r, t), (42) 

J 3 

where N c (t) is the number of particles in the conden- 
sate, and <pi(r,t) (i — 1,0, —I) is the condensate wave 
function, which is normalized as 

f dr£>i(r,t)| a = l. (43) 

i 

The condensate wave function is conventionally defined 
as <fii(r,t) = N c (t)ifii(r, t). However, throughout 
Sec. IIV1 the term "condensate wave function" refers to 
<Pi(*,t)- 

The field operator is then separated into the conden- 
sate and noncondensate parts: 

$ i (r) = <p i (r,t)a(t) + 6 i (r,t). (44) 

The noncondensate fraction is defined as 

_ N nc (t) _ N nc {t) _ J Y 



X ' N c (t) 
and, thus, we have 

8i(r,t) ln nc {t) 



tpi(r,i)a(t) 



n c (t) 



(fit(t)a(t)) 



(45) 



(46) 



We start with the Hcisenberg equation of motion for 
the operator cr (t)Si(r,t): 



ih 



dt 



(47) 



where T-L is given by Eq. ([2]). 

Expanding the right-hand side of Eq. (|4"7]) in powers of 
5i(r,t), and collecting terms of the same order of magni- 
tude, we obtain 

d f 4 

ih—(a*(t)$i(r,t))= ds^Rn(r,s,t), (48) 

n— 

where R n (n = 0, . . . , 4) is the sum of terms that contain 
the n-th power of Si (r, t) . The explicit expressions for R n 
(n = 0,1, and 2) are given below. (The terms R3 and 
i?4 are irrelevant when one makes an expansion up to the 
order of x 1 .) 

The first term on the right-hand side of Eq. (|47|) can 
be written in terms of the field operator ^(r, t) by using 
the following expressions: 



j 



— Qy(r,s,£) 



(49) 

^•(s,i). (50) 
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In the following, the argument t of ipi(r, t), <5j(r, i), a{t) is 
occasionally omitted to save space in lengthy expressions. 
By using Eq. ([2]), the last term in Eq. (j47| is rewritten 

as 



CO 



^^(8)^(8)^(8)^(8) 



3,1 



Cl 



e (/« wi («o$ oo^ ( g ) 



ot,j,k,g,l 



(51) 



This commutator can be calculated by using the relations 



&()",*), V>j(s,*)] = Q y (r,s,i), 



(52) 
(53) 



which follow from the orthogonality between the conden- 
sate and noncondensate. 



Substituting all the above results into Eq. (|47j) . the 
Heisenberg equation of motion can be rewritten as 



ihj t {a\t)Uv,t))= J dsj^ 

+ E 



d - - d 



(h )ji(pi(s)4>j(s)Si(r) + Q lj (r 7 s)a\h ) J ii)i(s) 
Co( - ^( s )'0]( s )V'/( s )V'j( s )4(r) + g ii (r,s)o t ^ t (s)| I (s)!/- J -(s)j 

<Ms)^j(s)^(s)^(s)<5 4 (r) + Qii(r,s)a t i>l(s)^i(s)i) k (s) 



Cl E (fa)jk(fa)gl 
a,j,k,g,l 



(54) 



Next, by substituting Eq. (|44|) into Eq. (|54j) and collecting terms according to the power of the noncondensate operator 
Si, we obtain i?„'s (?i = 0, 1, 2) in Eq. (|4"5)l as follows. 



- iflS ji-Q- t + (Mii + co(a f a - E IVfe(s) 



+ c i E (f<*)jk(f a ) g i(a'a- l)<p*(s)<pi(s)<pk(s) 



ih—tpjis,^ (s)a + £i(r) Q ij -(r,s,i)a t (5 i (s) ^ 

+ E 1 ~ ( ft -o)j;^(s)^(s)a t (5 i (r) + ^-(r, s^/io^O 3 ) 

+ c [ - </3*(s)^ z *(s)( ( 9i(s)(/3 : ,-(s)a t a t a(5i(r) + Qij(r, s)^*(s)ipj(s)a t a t a£ 7 -(s) 

+ Qij (r, s)v?/ (s)<^ (s)tfj (s)a t aa + (r, s)^,* (s)^- (s)a t a t a<5 ; (s) | 

+ c i E (f*)jk{fa)gl\ - ^(s)^*(s)</j ; (s)<y9 fe (s)a t a t a<5 i (r) + Qy(r,s)^ ; (s)(^ fc (s)^(s)a t aa 

a,j,k,g,l K 

+ Qij(r, s)v3*(s)(p fe (s)a t a t a(5/(s) + Qi,(r, s)< / 3*(s)<pi(s)a t a t a<5 fe (s) L 



(55) 



(56) 
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R 2 (r, B,t)=^2 ( s > *)<*j( s )<^( r ) 

3 

2tpi (s)tpi (s)<pj (s)<?j (s)a + a5i (r) - 99* (s)^* (s V; (s)a + o t ^ (s)& (r) 
Q ij (r,s)(^(s) l 5|(s)a t a l 5 / (s) + ^*(s)a t a t i5i(s)(5 i (s) + ipi(s)Sj{s)a' i aS j (s)] 



^2 I - {h )ji<pi(s)dt(s)5i(r) +co 



Cl (/a)jfc(/a) fl iS 
aj,k,g,l K 



2ip*(s)(pi(s)<f k (s)S](s)a ii a5 l (r) - V3*(s)^*(s)< ( 5 ; (s)a t a t (5 fe (s)(5 l (r) 



(r,s) (^*(s)a t a t <5;(s) l 5 fc (s) + (^ / (s)(Sj(s)a t a<5 fe (s) + ( / j fc (s) l 5j(s)a t a(5 ; (s) 



(57) 



To make a systematic expansion in powers of y/x, 
we expand the condensate wave function ipi(r,t) and 
the generalized noncondensate operator, defined by 
Ai(r,£) = tf(t)8 % {r,t)/yJW, as follows: 



^ 2> (r,t) 



Vi (r, t) = ^ 0) (r, t) + (r, t) +x } (r, t) + . . 

v „ ' 

¥>?>(«•,*) 



(58) 



Af (r,t) 



A, (r , i) = A 2 (0) (r, t) + ^ (r , t) + X 5kf ] (r, t) + 



Af (r,t) 



where we define 

tp\°\r,t) = ]im(pi(r,t), 

5<pV(r,t) = Km(^(r,t) - ^ (0) (r, i))/Vx, 
x->o 



(i) 



(0) 



(i) 



(59) 

(60) 
(61) 
(62) 



and so on. Note that the perturbative expansions in 
Eqs. (|58|) and (|59|) hold only if the system does not un- 
dergo a quantum phase transition as the total number 
density n (oc x 2 ) is increased from zero to the final value, 
during which the order parameter and energy spectrum 
change smoothly with ^fx- 

By expanding both sides of Eq. (|47|) up to the order of 
X°i X 1 ^ 2 an d X 1 successively, and using the orthogonality 
relation (see Appendix [A] for the derivation) 

(at(t)^(r,t))=0, (63) 

we obtain the equations that must be satisfied by the 
condensate wave functions at different orders ipf^ (r, t) , 
^ (r, t), i J 9 l - 2 ' ) (r, i) and the lowest-order noncondensate 

operator A^(r, i). Here, we outline the procedures 1-3 
for deriving those equations. 



1. Expansion up to the order of x°- 
By using Eq. (|63|) in the left-hand side of Eq. ([58)1 and 
neglecting the terms that contains j, (r, i) , we obtain 
(J dsi?o(r, s, t)) = 0, which leads to the time-dependent 

GP equation for ipf^ (r, t) (see Appendix [B] for the deriva- 
tion): 



E 



ihSj 



dt 



(h )ij + c NS i:j ^2 



(0)|2 



^ 0) (r,t) 



a,j,k,l 



(64) 



Here r/(°>(t) is an arbitrary real function relat- 
ing ipf*^ to ipf " 1 through a unitary transformation 



ipf Y (r, t) = ipf (r, i) exp i J dt'r/^ (t') /H . The dynam- 

L J 

ics of (p.- (r, t) is governed by the equation that is similar 
to Eq. (|6"4"|) but with the right-hand side being replaced 
by 0. 

2. Expansion up to the order of x 1 ^ 2 - 
Similarly, we can obtain the equation for the condensate 
wave function at the next order, y>^(r, t). It turns out 
that the equation that must be satisfied by </?j (r, t) is 
the same as that for tpf\r,t), i.e., Eq. (fM)) . In other 
words, the condensate wave function does not change 
at this order. Also, at this order, the equation of mo- 
tion for the noncondensate operator at the lowest order, 
A\°'(r,t), can be obtained by expanding both sides of 
Eq. (|48|) up to the order of x 1 / 2 . It is the time-dependent 
BdG equation (see Appendix [Cl for the derivation): 



ih— A 

dt 



(0) 



;(°)t, 



(T,t) = Y l [A ij A? ) (r,t)+B ij Ay»(r,t) , (65) 

3 
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where 



Aj =(ho)ij + 



(0)|2 



k.l 



AT (0) (0)* 



ciNj2(f a )kh(U) 



(o)* (o) 



fc,J I 



ft ,g 



i(0) 



(66) 



AT (°) (°) 



ciiV^(/ a ) fe ,(/ a ) is ^ V s 0) 

h,g 



(67) 



Here, is the projection operator onto the subspace 
orthogonal to the condensate wave function at the lowest 
order, (p^° (r). Its action on an arbitrary vector compo- 
nent fj(r) is given by 



Oi 0) °/i(r)=^|dr'Qg)(r,r')/ i (r') 



(68) 



where QgV.r') = %<5(r - r') - ^ 0) (r)^ 0) *(r'). For 
uniform systems, however, the elementary excitations are 
plane waves with nonzero momenta, which are orthogo- 
nal to the condensate wave function, and, therefore, the 
projection operator Q^y can be omitted. 

3. Expansion up to the order of x 1 - 
By using Eq. (|63[) and keeping the terms on the right- 
hand side of Eq. (|48| up to the order of x 1 , we obtain the 
generalized GP equation for the condensate wave func- 

(2) 

tion at this order ip\ (r), in which the effects of both the 
noncondensate number density and the anomalous aver- 
age are included (see Appendix [Dl for the derivation): 



E 



^% + (^o) 4J +coiV c %El^ 2) | 2 



CO 



Cl 



(2) . -* (2) . ~{R) (2)* 

n k iip) +n jk ip] +m}j '<p\ 



Here, A^ c < A^ is the number of atoms in the condensate, 



is the renormalized anomalous average, and .Fi(r) 



is defined as 



^i(r)=y ds|c iV c (j2\<P k °\s) 



E ( r > S )^ 0) ( S ) + ( r ' S ) Vf } * (s) 



CliV C E (/a), fc (/ Q ) S ^ s 0) *(s)^ (0) (s) 
a,j,k,g,l 



n*j (r, s)^ 0) (s) + (r, s)^ 0) * (s) 



(70) 



The matrix elements of the noncondensate number den- 
sity hij (r, s, t) and the anomalous average rhij (r, s, i) are 
defined in terms of the noncondensate operator at the 
lowest order A.\°\r,t) as 



^■(r,s,t)EE (A,( 0)t (r,t)A^ (s,t)), 



;(«>)/ 



(71) 
(72) 



From the time-dependent generalized GP equation for 
the condensate wave function ^j 2 '' (r, £) [Eq. (|69|) ] and the 
time-dependent BdG equation for the noncondensate op- 
erator A| ' (r, i) [Eq. (|65|) ] , both the dynamics and ther- 
mal equilibrium properties of a spinor Bose gas can be 
obtained. 

Using the Bogoliubov transformation, the time evo- 
lution of the noncondensate operator Kf 1 (r, t) can be 
expressed as 

UfWV U (A) M)J + A [u^(r,t))> 

where b\ and &^ are the creation and annihilation oper- 
ators of the excitation mode labeled by index A, and the 
coefficients u\ X \r,t) and v^ X \r,t) are given by 



v ( l X) (r,t) j 



v ( l X) (r) j 



(74) 



Here, u\ X \r) and u l - A ' ) (r) are the solutions of Eq. ([34|) 
with yly and By replaced by those in Eq. (|57)) . and 
is the energy of the excitation mode A. 



B. A uniform Bose gas in thermal equilibrium 



We apply the results in Sec. IIV Al to a uniform 87 Rb 
condensate in thermal equilibrium near absolute zero. 
In thermal equilibrium, the condensate wave function is 

(2) 

time- independent ip\ (r) , while the occupation numbers 
of excitation modes are given by the Bose-Einstein dis- 
tribution 



1 



[eW- M ]/(k B T) _l 



(75) 
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Here, the chemical potential u is taken to be the eigenen- 
ergy of the condensate wave function within an error of 
the order of 1/N, where N is the total number of par- 
ticles. The number density and anomalous average of 
noncondensed atoms are given by 



hij(r,s) = 



E 

M>0 



{<4 A) »^(s)/( £ W) 

+ v ( l X \r)v^*( S )[f(e^) + l]} 
m ij (r, S )= J2 {«f'*(r)f (s)/(eW) 

e( A )>0 



(76) 



For the uniform system under consideration, the con- 
densate wave function ipi is spatially uniform, while the 
excitation modes take the form of plane waves: 



^ A) (r) = ^ k V k - r , 



j 

w W(r) 



^ k) e ik - r , 



(77a) 
(77b) 



where k is the wave vector, and v is an additional index to 
distinguish excitation modes. The term J-i(r) in Eq ([70]) 
then has the following form: 



(78) 



k 

e>0 



£>0 



i.e., the nonzero contribution arises only from the exci- 
tation modes with zero momentum and positive energy, 
and it is vanishingly small in the thermodynamic limit. 

The set of coupled equations concerning the conden- 
sate and excitations in thermal equilibrium is then given 
as follows: 

1. The GP equation for the lowest-order condensate 
wave function ip^ : 



(qi 2 + c n) (p)°> + ci ^(F a )(/ a ) 



(0) 
3 



(o),„(o) (79) 



where is the lowest-order chemical potential, F a = 
n^(Pi*{fa)ij'fj ( a = x iVi z ) are the components of 

i,j 

the lowest-order spin density vector, and ipf^ 1 is normal- 
ized to unity: 



E 



(0),2 



1. 



(80) 



2. The BdG equation for the excitation modes at the 
lowest order: 



Aij (k) Bij 

-m -AUk) 



(l/,k)N 



(81) 



where 

AvQt) = (e k + qi 2 + c n) <% + W f (^ 0) )* 



ci 



E 



{Fa)(fa)ij 



k.l 



(82) 



= con^fVf + E c Mf a )ik(f a )jW k 0) <Pi 0) ■ 



(83) 



Here, e k = fi 2 k 2 /(2M) is the kinetic energy of a single 
particle with momentum ftk. 

3. The matrix elements of the noncondensate number 
density hij and the anomalous average expressed in 
terms of the excitation modes: 



E/(^{^ k) ^r k) / (o) (^ k) ) 

+ vf [/ (0) (^' k) ) + l] }, (84) 
E/||{^ k ^f k) / (0) (e^) 

/ (0) (e^) + l] }, (85) 



1 J 



where /(°)( e ) = l/{exp[(e - ^°))/(k B T)] - 1}. 

4. The generalized GP equation for the condensate 
wave function ip^ at the order of x 1 (Note that tp^' — 
ip^ as shown above Eq. ([55)) ): 



ci 



E 



(2) 



fc,i 



(86) 



where F° = n° y^p* [fo^ijtpf' i an ^ n " C anc ^ are 



(2) 



given by Eqs. (|28p and (|2T)]) . Here, is the renormal 



ized anomalous average which is described in Sec. IIV CI 
below, fi^ is the chemical potential at this order, and 

(2) 

the order parameter ip\ is normalized to unity: 

Ei^ 2) i 2 = 1 - ( 87 ) 



5. The number equation for the condensate and non- 
condensate number densities: 
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Note that the generalized GP equation (|86|) for the 

(2) 

wave function tp\ depends only on the lowest-order non- 
condensate operator A^(r) via hij = (A^°^A^) and 

fhij = (A| 0) A^ 0) ) [Eqs. dnj, (El)]- This is because the 
condensate and noncondensate operators are different in 
the order of magnitude, as shown in Eq. CHI) . 



or equivalently, 



C. Renormalized Anomalous Average 

The anomalous average term m^r, r'), which is de- 
fined in Eq. (|72[) . represents pair correlation of noncon- 
densed atoms, and can be expressed in terms of the ex- 
citation modes as in Eq. (|76p . However, the summa- 
tion over all excitation modes in Eq. (|T6[) would give 
an unphysical divergence. This divergence results from 
the replacement of the exact interaction by a contact in- 
teraction. This replacement amounts to assuming that 
all short-distance effects of the exact interaction can be 
encapsulated in one parameter: the s-wave scattering 
length. The effects of all higher-order terms, which rep- 
resent the multiple-scattering processes involving virtual 
states with high energies, are implicitly represented by the 
s-wave scattering length. Because fhij is first-order with 
respect to the interaction, taking into account the effect 
of pair correlation of noncondensed atoms on the conden- 
sate represented by CoTOy, C\fhij, would lead to a double 
counting of the terms that are second-order with respect 
to the interaction. This gives rise to the above divergence 
in the anomalous average term. 

To avoid this double counting, we need to go beyond 
the Born approximation and express the s-wave scatter- 
ing length a up to second-order with respect to the bare 
interaction. By applying the Lipmann-Schwingcr equa- 
tion (see [IB], P- 125) to low-energy collisions between 
two particles with a contact interaction, we obtain 



9 n A 2e° 



(89) 



k<k c 



k<k c k 



(90) 



where g is related to the s-wave scattering length by g = 
4irh 2 a/M, while g is the bare contact interaction. Here, 
e k = ^ 2 k 2 /(2M), £7 is the volume of the system, and k c 
is the cut-off of the momentum. 

For spin-1 atoms, there are two s-wave scattering 
lengths clq and for the total spin -Ftotai = and 2 
channels, respectively, and therefore, the corresponding 
coupling constants are given by 



9o =9o + 



So vl 
fi A-f 2e°' 



k<k c 



92 =92 



gj 



E O~0> 



(91) 



where go = 47r/i 2 ao/M and <?2 = ^f^a^jM. 

The spin-independent interaction cq and spin- 
dependent interaction c\ arc then given by 



co = 

Cl 



go + 2g 2 
3 ' 

h - go 



(92) 



By collecting all second-order terms with respect to 
the interaction, we obtain 



co E ™if + Cl E (fahiUhifhtfVt = co^ihijtp* +ci ^2 (U)ij{f a )kifhji<Pi 



"13 1-3 

3 



a,j,k,l 



(co - c )N c l^f I <Pi + ( c i - C 0^ C E (/«)«(/») klfkfWj 

3 I Oi,j,k,l 



E( ~ , ffo + 2 S2 1 (0) (0) \ * 

3 \ k<k c k / 



ki cirhji + 



k<k c k / 



Here, in obtaining the last equality, we have replaced iV c and tpi in the second line of Eq. (|93|) by N and (p\°\ re 
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spectively. This replacement causes an error of the order 
of x 3 / 2 , and therefore, is justified up to the order of x 1 ■ 
If the mean-field ground state is in the polar phase, 



i.e.. 
by 



- (J?) 

m 00 



(0,1, 0) T , the matrix element jTOqq 



m 00 + c n 



2e° 



= c n 



E 



is given 



(c n) 2 - (eg. + c n - e k ) 2 
1 

2^ 



(94) 



A(p/(p 



0.25- 



0.1- 



O- 1 



n cm 



10 7 



10 8 



10 9 



10 1 ' 



10 11 10 12 



where £k = \J e k( e k + 2c^n) > eg. Here in the first line 
of Eq. (|94[) we used the fact that — &o for 87 Rb. From 
Eq. ((94J), we find that mffi > 0. 



D. Ground-state phase diagram 



FIG. 5: (Color online) Density dependence of the rela- 
tive shift of the order parameter from the mean-field value: 



Atp/tp = max|</?^ — tp\"' |/ max \tpf' | (squares). The dot- 

ted line Atp/tp — 0.1 gives an estimate of the threshold below 
which the perturbative expansion gives quantitatively reliable 
results. 



By numerically solving the set of coupled equa- 
tions (|79 p -(|88 p we have calculated the ground-state or- 
der parameter (i.e., at T = 0) of a spin-1 ferromagnetic 
BEC up to the order of x 1 - The parameters of the sys- 
tem are the same as those given in Sec. MI Bl namely, 
those of 87 Rb atoms. Before discussing the result, let us 
first evaluate the threshold of the total number density 
"thros, beyond which the result obtained by the x 1//2 per- 
turbativc expansion deviates so greatly from the mean- 
field result that the perturbation method no longer gives 
quantitatively reliable results. Here, we define a measure 
of the validity of the perturbative expansion as 



Atp 



max \<p 

i,q 



(2) 



<P 



(0) 



max | (p. 



(0) 



(95) 



The perturbative expansion is valid if Atp/tp <C 1. The 
estimation of the value of Jtthres can be made in the 
following manner: the large ratio of co/|ci| ~ 200 
brings about a significant effect of noncondensed atoms 
on the spinor condensate; the condition for the effect 
of noncondensed atoms to be small is therefore given 
by Atp/tp < 0.1 or c n nc / (\a\n) < 0.01 (note that 
n oc \<p\ 2 ); using the expression for the noncondensate 
fraction n nc /n = 8(?ia 3 ) 1 / 2 /(3v / 7r), we obtain the condi- 
tion n < nthrcs = 10 10 cm" 3 . We have also solved the 
coupled set of the generalized GP and BdG equations 
for various values of n to calculate Atp/tp. The result is 
shown in Fig. [5l from which we find that the x 1 ^ 2 per- 
turbative expansion is valid for n < n-thrcs = 10 10 cm' 3 . 

Figure |6] shows the q-dcpendcncc of the longitudinal 
and transverse magnetizations of the condensate at ab- 
solute zero with n = 1 X 10 10 cm -3 . The mean-field 
result is superimposed for comparison. From Fig. [6j we 
find that the phase boundary between the BA and po- 
lar phases lies at q = % = 2.05|ci|n. The first-order 



self-consistent approximation discussed in Sec. IIIII with 
the same atomic density gives qh = 2.02|ci|n. These 

( R ) 

results show that the anomalous average further 
expands the region of the BA phase from the result of 
the first-order self-consistent approximation. This can 
be understood by considering the solution to the coupled 
equations (|T9" |) - (j55j) for q > 2\ci\n. The lowest-order con- 
densate wave function, which is the solution to Eq. (|T9"j) , 
is the same as that of the mean-field ground state, and 
it is given by tp^ = (0, 1,0) T for q > 2\c\\n. Since the 
atoms are condensed in the i = state, the matrices 
hij and ihij are dominated by the matrix elements noo 
and moo, respectively. The higher-order condensate wave 
function (f^ 2 \ which is the solution to Eq. (|86[) . can then 
be obtained as: 



(2) 



'(0,1,0) T (polar) if£>2 

/ V(2-Q/8 \ 
y /(2 + Q/4 (BA) if i < 2, 

ivy (2 -0/8/ 



(96) 



where £ ~ [q— 00(^00+^00" )]/ (l c l I 17 -)- The phase bound- 
ary between BA and polar phases is, therefore, given by 
£ = 2, or q = % ~ 2|ci \n + Co(noo + m oo )■ At absolute 
zero, noo and to ou "' ) arc both positive (see Sec. IIV C|) . 
Hence, the anomalous average further enhances the shift 
of the phase boundary toward the polar phase region. 
Note that the perturbative expansion breaks down in the 
critical parameter region 2|ci|n < q < % because in this 
region the system undergoes a quantum phase transition 
from the polar to the BA phase as the total number den- 
sity n is increased from zero to the final value [see below 
Eq. (p32"T) ]. However, the value of F±/n c shown in Fig. 
[5] (indicated by the double arrow) is found to be consis- 
tent with the expansion of the BA phase region at least 
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... r^/n c 




-10 12 3 



FIG. 6: (Color online) q dependence of longitudinal mag- 
netization per condensate particle F^/n c and that of trans- 
verse one F1/n c for a uniform ferromagnetic BEC at T — 0. 
The parameters are those of 87 Rb with atomic density n — 
1 x 10 10 cm~ 3 . The squares show the transverse magneti- 
zation numerically calculated by using the x 1 ^ 2 perturbative 
expansion, while the solid curve shows the mean-field result 
given by Fl/n c = s/l - (q/2\c 1 \n) 2 Q(2\c 1 \n - q)Q(q), where 
Q(q) is the unit-step function. The longitudinal magnetiza- 
tion, which is given by F^/n c = Q(—q) and shown by the 
dashed lines, is the same for the two approximations. The fer- 
romagnetic, BA, and polar phases are shaded in light, medium 
and dark grey, respectively. The phase boundary between the 
BA and polar phases lies at q = qb = 2.05|ci|n. The point 
indicated by the single arrow shows the value of F]_/n c at 
q/(|ci|n) = 2. In the mean-field approximation, F±/n c = 
at q/{\c\\n) = 2. The deviation of this point from zero in- 
dicates how much the BA phase expands from the mean- 
field result. The point indicated by the double arrow shows 
the value of F^_/n c that lies in the critical parameter region 
2|ci|n < q < qb- 

qualitatively. 

Figure [7] plots the value of % at the phase boundary 
between the BA and polar phases for various values of 
the total number density n. We find that the effect of 
quantum depletion on the spinor order parameter be- 
comes more significant for higher atomic densities, which 
in turn, leads to a greater expansion of the BA phase 
from the mean- field result. This trend in the shift of 
the phase boundary is clearly seen, eventhough the x 
perturbation method no longer gives quantitatively reli- 
able results for atomic density above rothres = 10 10 cm -3 . 
From these results, we conclude that the quantum deple- 
tion significantly alters the mean-field ground-state phase 
diagram of the spin-1 ferromagnetic BEC. In particular, 
when the atomic density is larger than 10 10 cm~ 3 , which 
is the case with usual experiments 0, > the system 
should be treated as a strongly interacting Bose gas. We 



c? b /(|cjn) 




FIG. 7: Position qb of the phase boundary between the BA 
and polar phases versus atomic density n. The dotted line 
shows the mean-field value qb = 2|ci|n. The values of qb for 
n = 10 11 and 10 12 cm , for which the perturbative expansion 
no longer gives quantitatively reliable results, are plotted only 
to show their rough estimates. 



shall examine this region in a future publication. 



V. CONCLUSIONS 

We have studied the interplay between the condensate 
and noncondensed atoms in a spin-1 87 Rb Bose gas in 
the presence of a quadratic Zecman effect. First, to in- 
vestigate the effect of thermal fluctuations on the con- 
densation and magnetism of the system, we have ap- 
plied the first-order self-consistent approximation and ob- 
tained the finite-temperature phase diagram. We find 
that the system can undergo a two-step phase transi- 
tion for a certain region of the quadratic Zeeman en- 
ergy: as the temperature decreases, the thermal gas first 
enters a nonmagnetized superfluid phase (polar phase), 
and then a superfluid phase with transverse magnetiza- 
tion (broken-axisymmetry phase). That condensation 
and spontaneous magnetization occur at different tem- 
peratures is characteristic of spinor condensates. Fur- 
thermore, via coupling to the magnetized condensate, 
spin coherence of noncondensed atoms in different mag- 
netic sublevels emerges, leading to magnetization of the 
noncondensate. By investigating the temperature depen- 
dence of magnetization of the noncondensate, we find 
that the magnetization of the condensate and that of the 
noncondensate are antiparallel to each other over a broad 
range of temperatures, except T < 0.01 To, where To is 
the transition temperature for a uniform ideal scalar Bose 
gas with the same atomic density. For T < 0.01 To, they 
become parallel to each other due to quantum depletion. 
This remarkable feature of the noncondensate at ultralow 
temperatures makes a distinction from high-temperature 
atomic thermal clouds. 

In contrast to scalar Bosc-Einstein condensates 
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(BECs), in spinor BECs the effect of a small fraction of 
noncondensed atoms on the system's magnetism cannot 
be ignored. This results from the fact that a large ratio 
of the spin-independent to spin-dependent interactions 
can significantly magnify the effect of a small number of 
noncondensed atoms. To examine the effect of quantum 
depletion and that of the anomalous average on the mag- 
netism of the system at absolute zero, we have applied the 
perturbative expansion in powers of x 1 ^ 2 , where x 1S the 
noncondensate fraction, to a 87 Rb spinor Bose gas. From 
the result, we have found that even a very small noncon- 
densate fraction can make a significant modification of 
the ground-state phase diagram from the mean-field re- 
sult. We have also found that when the atomic density 
exceeds a threshold nthros ~ 10 10 cm~ 3 , the deviation of 
the order parameter from the mean-field result is so large 
that the applied perturbation method can no longer give 
quantitatively reliable results. Therefore, a system with 
a higher density, which is usually the case with current 
experiments, should be treated as a strongly interacting 
spinor Bose gas. This is an interesting subject for future 
experimental and theoretical studies. However, to deal 
with this exciting regime, the system must be cooled be- 
low temperature T q( j, at which quantum depletion starts 
to dominates. Although the ratio T q d/T becomes larger 
as the atomic density increases, it is generally of the order 
of Xqd/To ~ 0.01, which presents a challenge in cooling 
techniques. 
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Appendix A: Derivation of Eq. (|63p 

The condensate operator a and noncondensate opera- 
tors Si (r) can be expressed in terms of the field operator 
ipi(v) and the condensate wave function <pi(r) as 

o= j dr^^*(r)^(r), (Al) 
k{v)=Y,Qij^) = /^dr'Q 4J (r,r')^(r'), (A2) 

3 J 3 

where Qy (r,r') = 5ij6(r — r') — <pi(r)<Pj(r') is the projec- 
tion operator onto the subspace orthogonal to the con- 
densate wave function ifi(r). From Eq. (|A2|) . we have 

/dr2>?(r$(r)=0, (A3) 

i 

i.e., the condensate and noncondensate are orthogonal 
to each other. The commutation relations between the 
condensate and noncondensate operators are given by 

[a, 0*1=1, (A4a) 

&(r), fit] =0, (A4b) 

[6 i (r),6}(r')]=Q ij {r,r'), (A4c) 

the others =0. (A4d) 

Using Eqs. ([AT]) and fA"2]) . we obtain 



(a + (^(r,t)) = My dr'^^(r',0^(r',t)j (j ds^ MM(M)) 
dr' J ds V>3 ( r '> Wit ( r - s > *) (*j ( s > *)) 



Py(s,r',t) 



■ j ds^Q«(r )S ,t) (j dr'^^.(s,r',t)^ (r',t)j 



N c J ds^Q«(r,s,i) W (s,i) 



(A5) 
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Hence, Eq. (|6"5)) is proved. 



order of magnitude of terms in R\ is, for example, 



Appendix B: Derivation of Eq. (|64[) 



By expanding both sides of Eq. (|48|) up to the order of 
X , and using Eq. we have 



ds(i? (r,s,t)) =0. 
From Eq. (|55|) for Rq(t, s, t), we obtain 



/ ds^Q y (r )S )|^ 

J ] Ik 

+ c N6 jk J2\d 0) (s)\ 



d 

lh5jk di + <yh °^ k 



Cl N ]T (/ Q ), fe (/ Q ) 3i ^ 0) *(s)^ 0) (s)^ 0) (s) 

a>,k,g,l 



(Bl) 



= 0. 

(B2) 



it at; 



a^iV^r)) 




Here, the number fluctuation operator is defined as 
AiV c = tfa- (a+a) = a+a- N c , (A7V C ) 2 = ((AiV c ) 2 ), 
and for a macroscopic number of particles in the conden- 
sate N c ^ N > 10 6 , the above term can be neglected up 
to the order of x 1 / 2 because l/\//V c <C 1. 

Consequently, only the first term (i? (r,s,i)) should 
be retained in Eq. (|C1[) up to this order, and thus the 
condensate wave function ip^ (r) must satisfy the same 

equation as Lpf^ (r, t), i.e., the time-dependent GP equa- 
tion: 



From the definition of the projection operator Qy(r,s), 
we arrive at the time-dependent GP equation: 



E 



-iMij— + (*»o) i3 - + coNSij \^k } \ 2 



cin (f«Mf a Wrv\\ 

a,j,k,l 



(B3) 



where r]^ (t) is an arbitrary real function corresponding 
to a unitary transformation as described below Eq. (1641) . 
Here, the number of particles in the condensate N c — 



E 



(1)|2 



a,j,k,l 



(C3) 



With the condition Hm<p^(r, t) = pi (t, t) and the 
x->o 

normalization condition, it can be shown that rf 1 '^) = 

T)(°\t) and i/j^ (r,t) = <p^(r,i), i.e., the condensate 
wave function does not change at this order. 



N(l - x) is replaced by the total number of N extj tne equation of motion for the noncondensate 



particles N with an error of the order of x\ which can 
be neglected up to the order of x° 



operator at the lowest order Aj (r, t) is obtained directly 
by expanding both sides of Eq. (|48j) up to the order of 

yl/2- 



Appendix C: Derivation of Eq. (|65|) 



First, by using Eq. (|63|) and expanding both sides of 
Eq. (|4"5)l up to the order of x 1 ^ 2 we obtain 



ds(<JZo(r,s,t)) + (i?i(r,s,t))) = 0. (CI) 



To calculate the second term in Eq. (|C1|) . we use the 
property of the condensate that the atomic number fluc- 
tuation in the condensate is of the order of AN C /N C ~ 
1 / \/N c . The expectation value (r, s, t)) then vanishes 
because compared with the lowest-order terms in Ro, the 



dsi? (r, s, t) 



'N c 
1 



/ N C 



(*i) 



-L= / daRi(r,s,t). (C4) 
'iV c J 



(*2) 



The first term on the right-hand side vanishes because 
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ip^\r,t) satisfies the ordinary GP equation, so 



The second term can be written as 



(*l)=-^=at a J ds^Q tJ (r,s)|^ 



{h ) jk + coN6 jk Y,\^\^ 2 



+ Cl N J2 (/«)j*(/«) ff «^ 1) *(8)^ (1) (s)^ 1) (8)i 

= 0. (C5) 



J 



(*2)=x; 



+ Cl 7V ]T (/a)y(/«Wi 0) >M>)Af (r) 



ds £ Q y (r, s) J c ^]T^^ 



ciN Y, (/aU(/ Q ) 5i [^ 0) *(s)^ 0) (s)A ( (0) ( S ) +^ 0) (s)^ 0) (s)Af t( 8 ) 



(C6) 



By separating terms containing (r, i) from those containing A^°^(r, t), we obtain the time-dependent BdG equation 
for the noncondensate operator A[°\r,t): 



ihjkf ] (r, t) = A^-Af (r, t) + i^-Af 1 ( 



;(°)t. 



(C7) 



where 



Cl iV^(/ Q ) y -(/ Q ) fc ^i 0) *(r)^ 0) (r) 



,k,l 



£ Off o co7V^ Vr* + Cl Ar^(/ a ) feff (/ a ) w ^V(°) oQ 



(o) 



fc.i 



.</■'< 



(0)* 



(C8) 



Appendix D: Derivation of Eq. (| (>!)[) 



By expanding both sides of Eq. (|48|) up to the order of x\ and using Eq. (|63|) . we have 



/" ds((i? (r,s,i)) + (i?i(r,s,i)) + (R 2 (r, s, t)>) = 0. 



(Dl) 
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The first term in Eq. (ID1|) is given up to the order of x 1 by 

ds(i? (r,s,t)) = f ds^Q. tJ (r,s)J ^ 

J 3 I k L v 

+ cJ((Ncf)-N c )5 jk J2\^ 2) (s)\ 2 



^ 2) (s) 



Cl 



(((7V C ) 2 )-7V C ) J] (/oU(/o) 9; ^ 2) *(s)^ (2) (s)^ 2) (s) 



(D2) 



a : k,g,l 



The number fluctuation in the condensate satisfies (AN C /N c ) 2 ~ l/iV c <C x 1 (f° r ^ ^ 10 6 ), and can be neglected up 
to this order. Equation (|D2j) then reduces to 



J ds(i? (r,s,i)) =N C j ds^Q y (r,s)| £ 



ih5 ik-Q t + {h )jk + c N c 6 jk m 



(2) 1 2 



a,k : g,l ) 

Using the fact that AN C /N C ~ l/\/N c -C x 1 ^ 2 ! the second term in Eq. (|D1[) vanishes up to this order: 



The third term in Eq. (|D1[) can be rewritten as 



,1/2 

/TV 5 



«x J 



+ Cl 7V c £ (/ a ) Jfc (/ a ) ffl ^.( B ) v ^ 



(D3) 



(D4) 



| ds(i? 2 (r,s,t)) =/v c | ds£Qy(r,s)jco£[n«(s) V f^ 

+ ci 51 {fa)jk{fa) g i n g i{s,s)ip k {s) + h* kg (s, s)ifi(s) + rh kl {s,s)ip g 2) * (s) I - Ji(r), (D5) 

a,k,g,l ) 

where ^(r) is defined as 

Fi{v) = j ds<L/v c r£ lv, (0) (s)| 2 J E [fi«(r,s)vf(s) + m ij (r,s)<pf ) *(s) 



(D6) 



By substituting the above expressions into Eq. (|D1|) , and 
using the definition of Qy(r, s) in Eq. (|A2j) . we obtain 

the generalized GP equation for tp\ (r) given in Eq. . 
Note that .Fi(r) does not change under the operation of 
Qy(r,s), i.e., 

|ds£Qy(r,s).F J (s)=.F i (r). (D7) 
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